library(ggplot2)

# Set working directory to project directory
setwd(wd)

load("data/dataforstage2_logit.Rda")

# Issue voting over time (not in paper)
# Identify five issues with highest scores
# Total IV
df <- dat[order(dat$iv_changes_gov_abs,decreasing=T),]
df$issueshort <- NA
df$issueshort[1:5] <- c("Strengthen WEF", "Nationalization","Cancel WEF","Strengthen WEF","Nationalization") 
pdf("absovertime_logit.pdf")
ggplot(df,aes(year,iv_changes_gov_abs)) + geom_point() + theme_bw() + geom_label(data=subset(df,!is.na(df$issueshort)),mapping=aes(year,iv_changes_gov_abs + c(1.4,1.4,1.4,1.4,-1.4),label=issueshort)) + xlab("Year") + ylab("Absolute IV") + xlim(1954,2010)
dev.off()

# Individual data points
pdf("ivpublicsupport_logit.pdf")
ggplot(dat,aes(gw_fvr_dk,iv_changes_gov)) + geom_point() + theme_bw() + geom_hline(yintercept=0,linetype="dashed") + xlab("Issue Voting") + ylab("Public Support (%)")
dev.off()

# Get descriptives for Table A2
sum(dat$iv_changes_gov_sig==1,na.rm=T)/sum(!is.na(dat$iv_changes_gov),na.rm=T)
sum(dat$iv_changes_gains_sig==1,na.rm=T)/sum(!is.na(dat$iv_changes_gov),na.rm=T)
sum(dat$iv_changes_losses_sig==1,na.rm=T)/sum(!is.na(dat$iv_changes_gov),na.rm=T)

summary(abs(dat$iv_changes_gov)[which(dat$iv_changes_gov_sig==1)])
summary(abs(dat$iv_changes_gains)[which(dat$iv_changes_gains_sig==1)])
summary(abs(dat$iv_changes_losses)[which(dat$iv_changes_losses_sig==1)])

dat$Policy <- as.factor(dat$policy)
levels(dat$Policy) <- c("No","Yes")
dat$Policy <- relevel(dat$Policy,ref="Yes")
pdf("ivpublicsupport_imp_logit.pdf")
ggplot(subset(dat,!is.na(dat$Policy)),aes(gw_fvr_dk,iv_changes_gov,colour=Policy)) + geom_point() + theme_bw() + geom_vline(xintercept=50,linetype="dashed") + geom_hline(yintercept=c(-1,0,1),linetype="dashed",colour=c("grey","black","grey")) + ylab("Issue Voting") + xlab("Public Support (%)") + scale_colour_grey()
dev.off()



# Make graph considering whether issue voting is significant or not (positive and negative)
implementation1 <- vector()
for(i in seq(20,80,20)){
  implementation1 <- c(implementation1, mean(dat$policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==0)]==1,na.rm=T))
}
implementation2 <- vector()
for(i in seq(20,80,20)){
  implementation2 <- c(implementation2, mean(dat$policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov>0)]==1,na.rm=T))
}
implementation3 <- vector()
for(i in seq(20,80,20)){
  implementation3 <- c(implementation3, mean(dat$policy[which(dat$gw_fvr_dk>i-20 & dat$gw_fvr_dk<=i & dat$iv_changes_gov_sig==1 & dat$iv_changes_gov<0)]==1,na.rm=T))
}

df <- data.frame(gw_fvr_dk=seq(20,80,20),imp=c(100*implementation1,100*implementation2,100*implementation3), iv=c(rep("None",4),rep("Opp-to-Govt",4),rep("Govt-to-Opp",4)))
df$gw_fvr_dk <- df$gw_fvr_dk-10
df$iv <- as.factor(df$iv)
df$iv <- factor(df$iv,levels=c("Opp-to-Govt","Govt-to-Opp","None"))
pdf("typeiv_imp_logit.pdf")
ggplot(df,aes(gw_fvr_dk,imp))  + geom_point(aes(colour=iv)) + geom_line(aes(colour=iv)) + theme_bw() + ylab("% Implemented") + xlab("Public Support (%)") + scale_colour_grey() + theme(legend.title=element_blank())
dev.off()

# Models. Table A3
mod1 <- glm(policy ~ iv_changes_gov_sig*gw_fvr_dk + totalchanges*gw_fvr_dk,data=dat,family="binomial")
mod2 <- glm(policy ~ iv_changes_gains_sig*gw_fvr_dk  + totalchanges*gw_fvr_dk,data=dat,family="binomial")
mod3 <- glm(policy ~ iv_changes_losses_sig*gw_fvr_dk  + totalchanges*gw_fvr_dk,data=dat,family="binomial")
library(apsrtable)
apsrtable(mod1,mod2,mod3,stars="default")

# Figure A2a
m <- summary(margins(mod1,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_gov_sig"),]
pdf("absivme_logit.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Total IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Figure A2b
m <- summary(margins(mod2,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_gains_sig"),]
pdf("ogivme_logit.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Opp-to-Govt IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()

# Figure A2c
m <- summary(margins(mod3,at=list(gw_fvr_dk=c(10,20,30,40,50,60,70))))
m <- m[which(m$factor=="iv_changes_losses_sig"),]
pdf("goivme_logit.pdf")
ggplot(m,aes(gw_fvr_dk,AME)) + geom_point() + geom_pointrange(aes(max=upper,min=lower)) + theme_bw() + xlab("Policy Support (%)") + ylab("AME Govt-to-Opp IV") + geom_rug(data=dat,mapping=aes(gw_fvr_dk,0))
dev.off()



